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Abstract 


An intense research effort over the last few years has produced several competing and 
apparently diverse methods for generating meshes. This paper reviews recent progress and 
emphasizes the central themes that we can expect to form a solid foundation for future 
developments in mesh generation. 

Introduction ^ 

Although long recognized as a major pacing item [24,49], mesh generation has only 
recently achieved the long sought aim of making possible the calculation of complete aircraft 
flow fields. Even now mesh generation for complex configurations is not a routine task and a 
period of maturation is still required before meshes can be quickly and efficiently created for 
any type of aircraft shape. We are, however, at an important juncture where techniques for 
tackling the problem are well understood, if not yet fully developed. It is therefore a suitable 
time to take stock and consider which approaches are likely to stay the course and provide the 
basis for reliable and efficient flow calculation methods. 

Early methods for calculating transonic flow over airfoils and other simple two 
dimensional shapes were based on conformal mapping techniques [2,3,20,43,70] and simple 
shearings [5,77]. At that time mesh generation was regarded as a subsidiary part of the flow 
algorithm. Indeed the problem was essentially that of finding a suitable coordinate 
transformation, expressing the flow equations in the curvilinear coordinate system and then 
devising a solution algorithm for the transformed flow equations. The paper of Thompson et 
al. [88] was a significant break from this tradition. It advocated the idea of using elliptic 
equation methods to generate meshes around arbitrary shapes. Although this approach had 
been tried before [104], Thompson et al. were the first to realize its full potential and to 
exploit the idea as a general technique. The paper had two other important effects. First, it 
focused attention on mesh generation as a problem in its own right, largely independent of the 
flow solver. Second, their work also recognized the need to move from a global description in 
terms of an overall mapping to a local viewpoint based on a mesh defined by a set of points 
and an ordering of the points corresponding to the coordinate directions. There is, of course, 
no need to use global mappings, and finite difference formulae can easily be constructed for 
numerically generated meshes. What is interesting, however, is that the apparently more 
elegant approach, of using coordinate transformations to achieve a global mapping of the flow 
equations, did not have sufficient generality to survive as a method for handling complex 
shapes. At the other extreme, the use of a non-aligned Cartesian mesh [19,65,103] once held 
the prospect of treating arbitrary geometries at the cost of finding adequate interpolation 
procedures for applying the solid wall boundaiy conditions. Unfortunately this cost has proved 
a formidable obstacle and the non-aligned mesh approach has won few adherents. The lesson. 
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which is often repeated, is that mesh generation eschews the most elegant solution without 
abandoning elegance altogether. In this case the preferred approach is a surface conforming 
mesh and a flow solver that uses only information about the mesh point positions. It is a 
compromise between the beauty of global coordinate mappings which prove too rigid and the 
excessive generality of a non-aligned mesh which is too difficult to couple with the flow 
solver. In this context it is important to note another early development, the introduction of 
the finite volume approach [45,66] that treats each mesh cell as a control volume to 
approximate the integral form of the flow equations. This is an inherently more accurate 
approach for discretizing the flow equations, a natural way of ensuring conservation and is 
narticularly well suited to the general surface conforming mesh. These concepts are now so 
widely accepted that it is easy to forget that there was a time when the choice was not so 
obvious. 

The paper of Thompson et al. [88] stimulated a great deal of work using elliptic 
equations to generate meshes in both two and three dimensions. The possibility of basing 
numerical mesh generation on hyperbolic and parabolic equation sets has also been explored. 
One unfortunate aspect of this emphasis on numerical mesh generation is that relatively little 
attention has been paid to alternative procedures. Algebraic methods [6,28,32,74], in 
particular, offer many attractive features. The work of Eiseman [28] is perhaps the best known 
example of this approach. He constructs a set of coordinate surfaces between two boundary 
surfaces and then defines an interpolatory function to provide a smooth variation in the third 
coordinate direction. The greater the number of intermediate surfaces used, the greater the 
degree of mesh smoothness that can be obtained. 

A particularly powerful algebraic method is based on transfinite interpolation. This 
approach to mesh generation was first introduced by Eriksson [31,32], who has developed the 
idea into an extremely effective technique. In common with all algebraic methods, this 
approach is very fast compared with numerical mesh generation. The most significant feature 
of transfinite interpolation, however, is the ability to exercise a high degree of control over the 
mesh point distribution, particularly, over the slope of the mesh lines which meet the boundary 
surfaces. 
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Fig. 1 Multiblock Variants (taken from Kutler, Ref. 49) 
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The conflicting requirements, of control over mesh quality on the one hand and 
flexibility to handle arbitrary shapes on the other, have pervaded mesh generation throughout 
its development. Indeed the essence of mesh generation is the need to maintain as much of 
both features without seriously degrading either. An important step towards securing a 
harmonious balance between these two constraints is the use of a multiblock strategy. This 
concept is to break the flow field into several smaller blocks (essentially an ultra-coarse mesh) 
and then generate separate meshes in each individual block. This concept was first formulated 
by Lee et al. [50] and has rapidly gained wide acceptance. There are several variants of this 
approach (see figure 1), depending on whether there is continuity of mesh lines at the block 
interfaces (composite), whether the individual blocks are entirely independent and merely 
overlap one another (overlapped), or block interfaces but no continuity of mesh lines 
(patched). 

At first sight it would seem that a blocked mesh structure is further step away from 
simplicity and elegance towards messy complexity. However the experience of many 
researchers who have embraced multiblock in spite of this drawback, provides convincing 
testimony to its role as an indispensable component of mesh generation for complicated 
configurations. Again one has been forced to relinquish a global viewpoint (the use of a single 
mesh) in favor of a local perspective (several blocks each able to accommodate part of the 
flow field and boundary surface). 

The discussion has so far centered on structured meshes composed of quadrilaterals in 
two dimensions or hexahedra in three dimensions. The mesh is structured in the sense that 

there is an implied set of coordinate directions within each block. In order to treat complete 
aircraft, this structure can still prove an impediment unless a large number of blocks are used. 
One is then faced with the problem of defining the blocks and their interfaces, which becomes 
increasingly difficult as the number of blocks increases. A radical alternative to a structured 
mesh is the use of triangles in two dimensions and tetrahedra in three dimensions. There is no 
longer any inherent regularity in the mesh and one has made a complete transition from the 
global to a local description. This characteristic gives unstructured meshes maximum 
flexibility in treating complex geometries while retaining a high degree of control over mesh 
point distribution. Unfortunately the extreme generality makes the generation of unstructured 
meshes a singularly difficult task. In two dimensions several triangle mesh generation 
schemes are available, but in three dimensions tetrahedral mesh generation schemes are scarce. 
The first report of a tetrahedral mesh around a complete aircraft is given by Bristeau et al. [18] 
This was a remarkable achievement, which has perhaps not received the recognition it 
deserves. No details of the mesh generation are given in the publications by this group, an 
unfortunate omission that has no doubt contributed to the lack of adequate recognition. 

There appear to be two main approaches to unstructured mesh generation which have 
proved successful in three dimensions. The first is based on the Delaunay triangulation and its 
dual geometric construct, the Voronoi diagram. This idea has been successfully exploited in 
two dimensions for a variety of applications. For aerodynamic calculations, use of the 
Delaunay triangulation was initiated by Weatherill [99,100] and developed in three dimensions 
by Baker [7,8,47,48]. The alternative approach is the moving front technique, successfully 
applied in two dimensions by Lo [51] and recently developed in three dimensions by Peraire et 
al. [63] and Lohner [52]. 

The distinguishing feature of unstructured meshes is their ability to handle three 
dimensional objects of arbitrary shape and complexity. This is accomplished without the need 
for introducing blocks, mesh controlling functions and other artifacts, which plague structured 
mesh generation and prevent its evolution to a fully automated procedure. The criticism that 
has been leveled at unstructured meshes is the difficulty of constructing efficient flow solvers. 
It is interesting to observe that the finite element method, which has exploited unstructured 
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meshes from its inception, has seen little development of the algorithms used to invert the 
mass matrix and solve the problem. Developments in Computational Fluid Dynamics began 
with simple structured meshes and this proved to be a stimulus for innovative ideas, which 
have led to highly efficient flow algorithms. Although it appears unlikely that implicit 
schemes can be adapted in any simple way to unstructured meshes, advances in explicit 
schemes (e.g. multi-stage Runge Kutta methods, local time stepping etc.) carry over equally 
well to triangular meshes. Indeed Jameson's unstructured flow solver [47,48] rivals the best of 
the structured flow algorithms for computational speed. Recent work by Mavriplis [55] 
indicates that multigrid will also work very effectively on unstructured meshes. 

It is likely that future developments will use combinations of both structured and 
unstructured meshes. In fact Nakahashi and Obayashi [59] have used such a combination, 
exploiting a structured mesh near solid boundaries to solve the Navier Stokes equations and a 
triangular or tetrahedral mesh in the remainder of the flow field. Actually the reverse 
procedure would be more appropriate for general geometries. It is the near field of closely 
coupled components such as wing/strut/nacelle and wing/store combinations that prove 
difficult when one tries to fit hexahedral cells. The likely outcome is that a tetrahedral mesh 
will be employed in such regions with an structured mesh covering the surrounding space. 

In the remainder of this paper we describe some of these ideas in more detail. Several 
excellent reviews and conference proceedings [30,73,85,86,89] have recently appeared and an 
AGARDograph containing contributions from several researchers will soon be published. 
From these publications the reader can get a comprehensive view of the current status of mesh 
generation. The text of Thompson et al. [90] is an excellent reference for much of the 

underlying theory. The aim of this paper is to provide a personal view by summarizing the 
most important developments and by attempting to interpret the many diverse ideas and 
techniques that have appeared. Also included is a brief discussion of two areas that can be 
expected to become increasingly important in the future. First, there is the question of mesh 
quality and, in particular, how changes in cell shape, aspect ratio and mesh stretchings can 
affect solution accuracy and algorithm performance [36,37,40,64,67,90,92,94,95,98]. These 
are difficult questions and our present knowledge is rudimentary. But it is likely that this area 
will receive much more attention, now that the major goal of complete aircraft meshes has 
been achieved. 

The second area is the use of solution adaptive meshes. This is also a field of growing 
interest and several papers on this subject have already appeared 

[1,14,17,25,39,41,42,54,57,58,60,61,63,83]. It is evident that the flow field around a 
complicated object such as a complete aircraft has a complicated structure and widely varying 
length scales. In order to resolve such features as shed vortices, shock wave patterns and 
ultimately separated flow, solution adaptive meshes will be needed no matter how many 
gigawords of computer memory become available. It is against this yardstick that further 
developments in mesh generation should be measured. How easily can the method cope 
automatically with an arbitrary shape, first generating an initial mesh and then adapting the 
mesh to the evolving solution so that all the salient flow features are crisply captured? 

2. Mesh Generation Methods 

In this section we provide thumbnail sketches of those mesh generation techniques that 
can be applied to complicated three dimensional shapes. For the purpose of illustration, 
formulae for the two dimensional implementation will be presented, and we therefore consider 
a transformation from physical coordinates x,y to another coordinate system defined by £,tj. A 
point in the flow field can then be associated with the position vector 

r(£/n) = (x(^,ti), y(^Ji)) 
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The reader should consult the references for a detailed derivation of these methods and their 
generalization to three dimensions. 

2.1 Numerical Techniques 

2.1.1 Elliptic Systems 

In two dimensions the Laplace system 

% +£ =0 

^xx yy 

( 1 ) 

TJ + T1 =0 

•xx 'yy 

may be solved to determine the coordinate transformation x(^,T|), y(^,rj). In practice the 
problem is inverted to solve for x and y as dependent variables [88,90]. Equations (1) then 
assume the form. 


°“« - P x 5ti + = 0 

«y^ - Py^, + Wtro - ° 


where 


2 2 
a = x^ + y^ 

P “ x 4 x r, + ^ 

2 2 
Y = x r| + y ii 

For a boundary conforming mesh, the boundaries are the lines % = ^ and tj = T]i, TI 2 say, 

and the solution of equations (2) thus defines a mapping from a rectangle in (^,q) space to the 
required region in physical space. The lines £ = const, and r] = const, correspond to mesh 
lines in physical space and the extremum principle for harmonic functions ensures that mesh 
lines do not overlap. 

2.1.2 Control Functions 


Unfortunately the Laplace system has a strong smoothing effect and often produces an 
undesirable mesh point distribution. This is particularly evident in regions of high curvature 
on the boundaries. Control of mesh point distribution can be accomplished by the introduction 
of controlling functions to generate a system of Poisson equations 


^xx + ^yy P ^ ,T ^ 

"Hxx + ^yy “ 


(3) 


The source terms P and Q are constructed to either attract or repel points about a certain 
position or line. For example 

P(^,rj) = -a sgn(^ 0 )e _c I I , a > 0, c > 0 (4a) 

attracts lines towards the line while 
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(4b) 


P&Tl) = -a sgn(^ 0 )e^ c[( ^ o)2 + ^ 

attracts t ,— lines to the point (^o.'Ho)- Ways of automating the choice of the amplitude a and 
decay factor c have been suggested, and more general forms of the source terms are available 
[56,79,90]. However, non-overlapping of mesh lines can no longer be guaranteed and the 
specification of control functions requires some care. 

2.1.3 Hyperbolic Systems 

Other sets of partial differential equations might be expected to produce satisfactory 
meshes. An interesting development of Steger and Sorenson [80] along these lines is the 
solution of the following equations 


x $ x q + y l; y q =0 
x £ y q “ x q y £ = V( ^ 


(5) 


which can be shown to define a hyperbolic system. The first equation of this pair is the 
orthogonality condition; the second determines the local cell area according to a specified area 
distribution V(J;,q). If the x and y coordinates of one boundary surface are defined on q=qj 
say, then the system can be marched out in the direction of increasing q. The outer boundary 

surface cannot, of course, be specified for a hyperbolic system. This is often quite acceptable 
for external flow problems, but makes a hyperbolic system unsuitable for generating meshes 
for some geometries, typically those used in internal flow problems. 


2.2 Algebraic Methods 


2.2.1 Multi-Surface Fitting 

An interesting approach that was introduced by Eiseman [28] is based on the idea of 
interpolating curves along specified directions rather than through specific points. Consider 
two boundary surfaces ri(|) and r N (£) corresponding to the lines T|=r|i and q=q N - Let 

r(^,q) = (x(£,q), y(^,r|)) be the position vector of the point (x,y) and introduce intermediate 
surfaces r* (q), i=2,...,N— 1 which are constructed to direct the q— lines from the boundary 
surface q=qi to the opposite boundary surface q=q2- We now construct a curve that is tangent 
to the direction q j(£) - Ij(£) at the position q=q - t Thus we require 


3r 

3q = 


N— 1 


l A iV i <’lXi i+1 ©-q©) 

i=i 


( 6 ) 


where Ajare suitable normalizing constants and the \gj(q) are univariate interpolating functions 
subject to the constraints 


ViOlj) = §jj i,j=l,2,...,N— l 
On integrating eqn (6) we obtain 


N— 1 

r(5,'l) = i 1 ©+ l Aj 0 ,di) (ij +1 © -q©) 

i =1 


where 


GjOD 


r'n 

V:(s)ds . 

\ 


The choice rjf^) for the integration constant ensures that r(^,Tj) coincides with r^) when 
ti=tii. We can now choose Aj = 1/G i (r|jP to ensure that r(^,r() coincides with the outer 

boundary surface r N (0 at 

This approach offers great generality and a high degree of control over mesh point 
distribution. From a practical viewpoint it must be difficult to specify the intermediate 
controlling surfaces for regions contained by highly complicated boundaries (e.g. the space 
around a nacelle/strut/wing). On the other hand, once one has accepted the need for a blocked 
mesh, a more specific interpolation procedure such as transfinite interpolation is probably 
adequate. 

2.2.2 Transfinite Interpolation 


The theory of multivariate interpolation, developed in general terms by Gordon [38], 

has been successfully exploited by Eriksson [31,32] as a means of generating meshes 
exhibiting a high degree of control over the mesh point distribution. Unlike Eiseman's 
multisurface approach, only information about the distribution of boundary points and the 
slope of mesh lines meeting the boundaries is needed. However, a careful choice of the 
interpolating functions is required to obtain a smooth mesh when one or more of the boundary 
surfaces contains a slope discontinuity. This aspect is discussed further in reference 31 which 
provides a clear description of Eriksson's work. A particularly attractive aspect of this 
approach is the way in which an interpolating function in two or three dimensions can easily 
be constructed from univariate interpolants. To illustrate this in two dimensions we again 
consider a mapping of the rectangular region in ^,rj space defined by the boundaries £=^, 
^=^ 2 . T|=r|i and T|=T|2- We require the vector function r(£,r|) = (x(£,r|), y(£,r|)) given a 
prescribed functional variation r(£i,ri), rC^Jl), Th<ri<ri 2 and r(£,Tli), r(q,ri 2 ), £,i <^<^2 on the 
boundaries. First, we define two projections that correspond to a pair of univariate 
interpolation schemes, 


where 


E = I^pTl) 4^(5) + i(^ 2 dl) 4> 2 (S) 

\ E = KS/np YjOi) + r(^,T| 2 ) V 2 (r|) 

4*2^9 p = 0 » $2^2^ = ^ 

¥j(Tip = l, ¥Pti 2 ) = 0 

Y 2 ( t lp = 0 ’ V 2 0l 2 ) = l 


(7) 


( 8 ) 
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We now define the transfinite interpolating function as the Boolean sum of these projections 

r = (ic^ + - jc^) r (9) 

= r(^pTi) + r(^ 2 ,ri) <> 2 (^) + r^Jip Vj(r l) + r(^,r| 2 ) y 2 (n) 

-r^j.rip <t> 1 (^) Vj(ti) -rC^.Tip $ 2 ® VjC^l) 

-r(^ 1 ,'n 2 ) (J)^) V 2 (ti) - r(^ 2 ,ri 2 ) <t> 2 (^) y 2 (il) 

Along the boundaries the interpolating function reduces to the prescribed form. The functions 
<t>i, <{> 2 , Yi and V 2 can be chosen arbitrarily provided they satisfy the conditions (8). The 
simplest possible choice corresponding to linear interpolation in £, and r\ is given by, 

♦ 1 ® = • * 2 ® = If^ti 

The projections (7) can be generalized to allow for a specification of the derivative of r 
normal to the boundaries. The univariate interpolants <|> and \\f must then be osculatory 
functions that satisfy conditions (8) together with a further set of constraints on their 
derivatives at the end points. 

2.2.3 Sequential Mapping 

Another algebraic method that retains something of the spirit of the original coordinate 

transformation methods based on conformal mappings and shearings is the use of a sequence 
of mappings [6,44]. This approach has the virtue of retaining a high degree of control over 
mesh point distribution. By applying several simple mappings it is possible to reduce a 
relatively complicated configuration to a simple generic shape. 

For example, a wing/fuselage/tail combination can be mapped into a pair of surfaces 
plus part of the symmetry plane. A Joukowski mapping followed by a shearing will take the 
fuselage into the symmetry plane. Further combinations of conformal mappings and shearings 
can be constructed to reduce the tail to a single sheet and unwrap the wing. It is then fairly 
straightforward to interpolate a set of coordinate surfaces which conform with the mapped 
wing surface and contain the tail sheet. In addition the coordinate surfaces are required to 
coincide with the mapped fuselage crown line on the symmetry plane. The set of mesh points 
in mapped space is now passed through the inverse of the original mapping sequence. This 
restores the aircraft geometry and produces a mesh that conforms with all boundary surfaces. 
A very similar approach was pursued independently by Shmilovich and Caughey [72]. 

The technique is fairly powerful, and the surface mesh for two examples of meshes 
generated by Baker's method are presented in figure 2. This mesh generator was linked to 
Jameson's finite volume scheme to produce the first published Euler flow solutions over 
wing/fuselage/tail combinations [44,46]. Although this concept could, in principle, be 
extended to include engine nacelles, it becomes increasingly difficult to fit the coordinate 
surfaces in mapped space, and some degree of user intervention would then be required to treat 
different geometries. As we mentioned earlier, the use of global coordinate transformations is 
not sufficiently flexible to deal with complete aircraft configurations, leading one to accept 
that a multiblock approach is unavoidable. 
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(a) Boeing 747—200 


(b) Lear Jet 


Fig. 2. Surface Mesh from Baker, Ref. 6 and 44 


2.2.4 Grid Blending 

Recently Steinhoff [82] introduced the idea of blending a collection of meshes, each 
one of which is generated for a separate region, to form a smooth global mesh. In some 
respects this idea has some of the features of a multiblock method in so far as it is based on 
defining several relatively simple, separate meshes rather than a global transformation. 
Suppose, for example, that there are N meshes and that the vector 


t,jk< m) = <V m >’ yij k <m >- V" 1 ” 


represents the point with indices (i,j,k) generated by the mth mesh where m=l,2,...,N. The set 
of points represented by 1 ^( 1 ) might be the mesh around the wing, 1^(2) the mesh around 

the fuselage, ^^(3) the mesh around the tail and so on. The actual point r^ corresponding to 
the indices (i,j,k) is then formed by the weighted combination 


N 


I V m > V m) 


f _«i=l 

hik~~w~ 


( 10 ) 


i 


m=l 


V m) 
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dominates in the region where it is needed. The weights for the other meshes should be small 
in this region, and decay to zero near the boundary surface where the mesh is determined 
entirely by the component mesh r^Cm). 

Although this approach is attractive, the choice of the weighting functions is by no 
means easy and will presumably have to be done on a trial and error basis for each new 
configuration. How difficult or easy this task is compared with the construction of a blocked 
mesh will determine whether the blending mesh approach is a viable technique for treating 
arbitrary complex configurations. 

3. Multiblock 

The need to handle complicated geometries while maintaining adequate control over 
mesh point distribution and cell shape can perhaps be best achieved by introducing a 
multiblock structure. The flow field is broken up into a number of blocks by defining a set of 
surfaces which will represent the block boundaries. The individual blocks can then be meshed 
and the flow algorithm constructed to exploit the blocked structure. In fact the flow solution 
can proceed in each block independently, with information between contiguous blocks being 
passed at the common interface of two block boundaries. This makes multiblock particularly 
suitable for computer architectures with a parallel processing capability. Various possibilities 
arise depending on what degree of continuity is required at the block interfaces (see figure 1). 

3.1 Overlapp ed 

If we do not define a precise interface but introduce separate meshes for each 
component, we obtain a system of overlapped meshes. This approach was originally 
considered by Atta [4] and has been extensively developed by Benek et al [10—12]. The lack 
of any constraint at the block boundaries makes mesh generation for the individual blocks 
much easier. However the penalty for this facility is the problem of transferring information 
from each component mesh to its neighbor. This difficulty has some similarity with the 
problem of applying solid wall boundary conditions on a non-aligned mesh. Similar 
reservations about accuracy of the interpolation procedures and the ease of maintaining 
conservation apply to overlapped meshes. Recent work by Berger [13] shows that 
conservation can usually be achieved though with some difficulty, particularly in three 
dimensions, and it is by no means clear whether stability could be assured in the general case. 

3.2 Patched 

If we define a block structure and require the separate meshes to conform with the 
surfaces of their respective block boundaries, we obtain a patched mesh. There will, in 
general, be no continuity of mesh lines from neighboring blocks at the block interfaces. 
However, interpolation at the interfaces is now less demanding than that required by an 
overlapped system, and this approach retains the advantage of allowing a highly refined mesh 
in specific regions without imposing unnecessary refinement elsewhere. A patched mesh was 
used by Baker et al. [9] to obtain a refined mesh in the vicinity of the tail for a 
wing/fuselage/tail combination. This technique has also been successfully applied to solve the 
Navier-Stokes equations around an F-16 [34,75,76] (see figure 3). 
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Fig. 3. Mesh from Sorenson, Ref. 75 


3.3 Composite 

The composite method can be regarded as a special case of the patched approach in 
which mesh lines are required to be continuous across block interfaces This has the drawback 
that mesh refinement in one block, involving an increase in the number of mesh points on 
block boundaries will induce a corresponding refinement in neighboring blocks and so on 
throughout the entire mesh. However, it is probable that the advantage of mesh line continuity 
at the block interfaces outweighs the disadvantage of requiring an unduly fine mesh in certain 
regions. Mesh smoothness is further enhanced by requiring slope continuity as well, and the 

extra burden this places on the mesh generation within each block is amply justified by 
improved accuracy that results. 



Fig. 4. Surface Mesh 
from Yu et al., Ref. 106 
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3.4 Miscellaneous Results 


Over the last two years several examples have appeared in the literature showing flow 
calculations on structured meshes for aircraft geometries. Yu et al. [105,106] have used a two 
block Poisson solver to obtain a C— H mesh around a high tail/aft mounted propfan 
configuration (figure 4). Eberle and Schwarz [27] have generated a single block H— H mesh 
around a wing/canard/fuselage/tail combination (figure 5). They used an elliptic mesh 
generator that solves the biharmonic equation, thus gaining a higher degree of smoothness than 
would be obtained from the Laplace system. A similar configuration, but without a vertical 
tail, has been treated by Eriksson et al. [33], who used transfinite interpolation to produce a 
two block H-O mesh (figure 6). 




Fig. 6. Mesh from Eriksson et al., Ref. 33 
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A mesh generator that is based on only one or two blocks must necessarily be restricted 
to configurations of one generic type. Other types of configuration could be handled only after 
extensive recoding of the mesh generation software. More general methods based on several 
blocks have been developed by Weatherill and Forsey [101] and Shaw et al. [71], who use a 
Poisson solver in each block (figure 7), and also by Fritz et al. [35] and Seibert [69] (figure 8), 
who generate an initial mesh algebraically and then iterate with a Poisson solver to obtain a 
smooth point distribution. 



Fig. 7. Mesh from Weatherill and Forsey, Ref. 101 
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The inclusion of nacelles, struts and stores introduces an added degree of difficulty. 
The work of Vigneron et al. [96,97] (figure 9), who used algebraic generation followed by a 
Poisson type smoothing is therefore particularly impressive. Another outstanding contribution 
is the paper by Sawada and Takanashi [68]. They used transfinite interpolation to produce a 
mesh around a complete aircraft with over— wing nacelles (figure 10). 



Fig. 9. Mesh from Vigneron et al., Ref. 97 


3.5 Further Developments 

Once one has accepted the need for a multiblock scheme and the extra data structure 
that this entails, it makes sense to keep the individual blocks fairly simple (i.e. as close to 
rectangular as possible) in order to simplify the mesh generation process within each block. 
To treat a complete aircraft configuration may require around 100 blocks, but this can certainly 
be accommodated if the data structure is sufficiently general and contains complete 
information about the block face and edge contiguities and the relative orientation of 
neighboring blocks [15,16,81,87,101,102]. Transfinite interpolation, with prescribed mesh line 
slopes at the block boundaries would appear to be a perfectly adequate method of mesh 
generation within each block. This could be followed by a few iterations of a Laplace solver 
if further smoothness were needed. The remaining problem is the decomposition of the 
flowfield to give the block structure. Ideally this should be completely automated, but this 
may be difficult to achieve for arbitrary aircraft shapes. The consequent need for good 
interactive graphics has been mentioned by several authors and it appears that structured mesh 
generation for highly complex geometries will continue to rely on the intervention of a skilled 
user. 
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Fig. 10. Mesh from Sawada 
and Takanashi, Ref. 68 


4. Triangulations 

A review of unstructured mesh generators prior to 1980 has been given by Thacker 
[84]. Significant developments, however, have taken place since then and efficient 
triangulation methods are now available for two and three dimensional problems. It is worth 
noting that any planar triangulation possesses a degree of structure that is not present in three 
dimensions. First we observe that given a set of V points in a plane such that B lie on the 
convex hull, the number of edges E and the number of triangles (or faces) F are specified by 

the formulae 


E = 3(V - 1) - B 
F = 2(V — 1) — B 


In particular, since B = O(V^) 


we see that for a large number of points, 


E ~ 3 V and F ~ 2V 
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Moreover, if k denotes the average number of edges meeting at a point then, since each edge is 
associated with exactly two points, the number of edges is also given by 


E = 


kV 

2 “ 


Combining this with the above expression for the number of edges leads to the following 
expression for k. 


k = 6 — 


6 2B 
V _ V" 


and so for large V we have 

k ~ 6. 


The average number of edges meeting at a point is therefore six, and since at least 
three edges must meet at any point, this places a fairly rigid constraint on the expected 
variation in the number of edges meeting at any point. In other words, we do not expect too 
much variation between triangulations for a given set of points in two dimensions. In 
particular, a high degree of regularity, in terms of the number of edges meeting at any point, is 
to be expected no matter how we triangulate the points. 

There is no such invariance in the number of tetrahedra, triangular faces and edges in 
three dimensions. It is well known, for example, that a cube can be cut into either five or six 
tetrahedra. In fact, it is possible to find triangulations of N points in three space which contain 
0(N) points and also triangulations which contain 0(N 2 ). Intuitively one would expect that an 
0(N) triangulation would be likely to contain better shaped tetrahedra and lead to a greater 
degree of mesh regularity. The important requirement is, therefore, to find a triangulation 
scheme and mesh point distribution that will achieve this aim. 


4.1 Delaunay Triangulation 


The Delaunay triangulation of a set of points and the dual geometric construct, the 
Voronoi diagram, are extremely fertile concepts that have been the subject of considerable 
theoretical investigation and have found numerous practical applications. The Voronoi 
diagram marks off the region of space that lies closer to each point than the other points. This 
is illustrated for the planar case in figure 11. The solid lines make up the Voronoi diagram 
which form a tessellation of the space surrounding the points. Each Voronoi tile (e.g. the 
hatched area around point P) consists of the region of the plane that is closer to that point than 
any other. The edges of the Voronoi diagram are formed from the perpendicular bisectors of 
the lines connecting neighboring points. In general three edges will meet at a vertex which 
must be equidistant from three forming points (e.g. points P, Q3, Q4 in figure 11) and hence 
each vertex is the circumcenter of the triangle formed by three points. This determines a 
unique triangulation known as the Delaunay triangulation and is such that the circumcircle 
through each triangle contains no points other than its forming points. 
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The complete mesh is generated by the triangulation of a cloud of points surrounding 
the aircraft. The point distribution can be defined in any way whatsoever, and it is this 
generality which makes the Delaunay approach so powerful. The surface geometry for the 
F-15 was based on a network of points originally generated for use by a panel method. 
Flow field points were introduced in three sets corresponding to a farfield, midfield and 
nearfield distribution. The farfield points were generated by defining a Cartesian box with 
a relatively coarse distribution covering the entire field. A finer distribution of points, 
again in a regular array, was added in a region extending a few wing chords around the 
aircraft. Finally, the nearfield points were generated by placing points along normals 
directed outward from the surface points. Any points which fall inside the aircraft structure 
are detected and removed. The remaining points are triangulated to create a mesh which 
varies from a fine distribution of cells near the aircraft to a coarser distribution in the 
farfield. 

The two photographic plates show the surface mesh and a flow solution for the F-15. 
The mesh was generated for one half of the aircraft; the mesh shown in plate 1 and the 
surface contours shown in plate 2 have been reflected about the aircraft plane of symmetry. 
The mesh surrounds one half of the full configuration, including the interior of the duct 
which extends from the fuselage inlet to the nozzle at the rear. Of the 77,000 points in the 
mesh, approximately 5,000 are on the aircraft surface, and the complete set of points has 
been connected to form a collection of 460,000 cells. The surface triangulation shown in 
plate 1 provides a good indication of the mesh resolution, and plate 2 shows the pressure 
distribution for a case run at a freestream Mach number of 0.81 and an angle of attack 
of 4.84 degrees. In the color coded picture (plate 2), red represents low pressure and blue 
represents high pressure with the fringe bar indicating the graduations in between these 
extremes. 
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Plate 1 
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Plate 2 
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These concepts generalize to higher dimensions. In particular, the Delaunay 
triangulation of three space is the unique triangulation such that the circumsphere through each 

tetrahedron contains no points other than its forming points. In two dimensions this circle 
criterion can be shown to be equivalent to the equiangular property that selects the 
triangulation which maximizes the minimum of the six angles in any pair of two triangles 
which make up a convex quadrilateral. No equivalent characterization is known in three 
dimensions but the circle criterion can still be regarded as selecting a good triangulation for 
the given set of points. 

The Delaunay triangulation has been applied to mesh generation for structural problems 
by Cavendish et al. [21] and for semi-conductor device simulation by Cendes [22,23] et al. A 
three dimensional, Delaunay based, mesh generator developed by Baker [7] has been linked to 
the finite element flow solver of Jameson and reference 47 contains the first published 
example of an Euler flow calculation over a complete aircraft. Detail of the surface mesh 
around the inner nacelle/strut/wing area of the Boeing 747—200 is presented in figure 12(a) 
and figure 12(b) shows the surface triangulation associated with a tetrahedral mesh that was 
generated for the F— 15. The generality and flexibility of this method is clearly evident and the 
technique is capable of generating tetrahedral meshes around any three dimensional object. 



(b) McDonnell Douglas F— 15 


Fig. 12. Surface Triangulation from Baker, Ref. 7 and 47 
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4.2 Moving Front 
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An alternative approach to the triangulation problem is the moving front technique first 
suggested by Lo [51]. This idea is illustrated in figure 13, where a layer of triangles has 
already been placed around an internal boundary surface. The outer edges of this layer forms 
the front and a new triangle is constructed on each edge by connecting the edge to the nearest 
point in such a way that the shape of the new triangle is acceptable. When triangles have 
been formed on all front edges, a second layer of triangles will be in place. The procedure is 
repeated with the new front until the whole set of points has been triangulated. 



The extension of this idea to three dimensions is not an easy task since the front, 
formed by a layer of tetrahedra, is now made up of a surface of triangular faces. It is 
presumably necessary to take care when determining the intersections of planar faces, and to 
ensure that no overlapping of tetrahedra occurs or that holes are left when the front folds over 
on itself. Significant progress has been made in this area, however, by Peraire et al. [62,63] 
and also by Lohner [52]. Figure 14 shows detail of the surface triangulation associated with a 
tetrahedral mesh that Peraire generated around a wing/canard/fuselage/tail combination. 



14. Surface Triangulation for Wing/Canard/Fuselage/Tail 
Configuration (taken from Peraire et al.. Ref. 63) 
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5. Mesh Type and Solution Accuracy 

An inevitable consequence of creating a surface conforming mesh is the introduction of 
mesh stretching. Some degree of mesh stretching is clearly desirable in order to achieve an 
efficient distribution of mesh points throughout the flow field. However, the stretching should 
be smooth to prevent spurious, numerically generated effects from contaminating the flow 
solution. Investigation of this subject is still in its infancy and an improved understanding of 
this important area is clearly necessary. Two aspects of mesh distortion are considered below. 

5.1 Truncation Error 

To illustrate ideas we consider a one dimensional stretching from the interval 0 £ £ < 1 
to physical space defined by the interval x Q < x < x^. In most cases there will not be an 

explicit mapping from % to x. The variation in mesh width can, however, be regarded as the 
result of an implied transformation x(£). We assume that the mapped interval is divided into 
N equally spaced increments A£, = 1/N and write xj = x(jA£) with the further requirement 

x o <x l <x 2 < - <x j < - <x N 


We also write hj = Xj — Xj_j 

and let h = max h; 
j 

As Thompson et al. [90] point out, it is necessary to distinguish between two senses of order. 
First there is the behavior of the error as the number of points in the field is increased while 
maintaining the same relative point distribution over the field. Consider the approximation to 
the first derivative f x given by 


f. 


j+1 ~ f j-l - + A^ 2 


•j+1 X j-1 6xe 






^||f^) + 0(A^ 4 ) 




(id 


The first term is, of course, the required derivative f x ; the second term is the principal part of 
the truncation error, which is formally second order. An alternative definition of order is 
based on the behavior of the error as the relative point distribution is changed so as to reduce 
the spacing locally with a fixed number of points in the field. Expanding Vi and fj_j as 


functions of x about the position Xj we obtain 


Vl 


in. 


x j+i 


— x 


= f 


j-1 


+ j( h ; 


‘j+l 


h j) f xx 


, lr.2 

+ 5 [h j+l 


- h j + l h j +h j ]f xxx + H0T 


( 12 ) 


But 


h j+l ” h j = X ^ A ^ + 0(A ^ 4) 

hj + l ~ h j+ ihj + hj = x^ 2 + 0(A^ 4 ) 
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Hence 


X^j - Xj_| - f x + (x ^ f xx + / f xxx ) + H0T 


(13) 


Now 


h. = A£ + O(A^) 


so we have 


!i±! ~ 


h? x t 


K^ =f ‘ + ^ ( f f ~ + * f - )+H0T 


(14) 


It is apparent that a necessary condition for the approximation to be second order in the second 
sense is that the term 


^ = 0(1) (15) 

This is discussed by Thompson et al. [90], who point out that a similar condition is required to 
maintain local order in two dimensions. A further condition must also be imposed to limit the 
rate at which the Jacobian approaches zero. This extra condition is essentially a limit on the 
degree of nonorthogonality that can be tolerated. 

It should be noted in passing that expression (13) or (14) is exact for linear f. This is 
not the case if we use the exact mapping derivative x^ in place of Xj + j — Xj_j. For then, we 

have 


Hit 1 






Etf + 
V, XX 


$ 


f ) + HOT 


XX 


It is therefore preferable to evaluate the metric coefficients numerically by the same difference 
representation as is used for the dependent variable. This conclusion was probably first stated 
by Steger [78]. 


The issue of local accuracy and stability on nonuniform meshes has been examined in 

h i+l 

some detail by Turkel [92]. He first defines a local stretching factor r. = f ' and then 

J J 

classifies three groups of stretchings: 

(a) Quasi-uniform (or algebraic) if 

r. = 1 + 0(h p ), p > 0 

(b) Exponential if (a) is not valid but 


1 + 0(h) and 


itu 

r j 

(c) faster than exponential if neither (a) nor (b) is true. 
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Now hj + j — hj = hj(ij — 1) and it follows that if the stretching is quasi— uniform then from 

expression (12), the approximation remains locally second order (i.e. second order in the 
second sense described above). A quasi— uniform stretching must therefore satisfy condition 
(15). Further we can write 


At 2 

h;,i + -i- Xte 

r . = 4±i = _£ ' — + HOT 

] h. A e2 

J 


= 1 + At + HOT 


or 


r. 

J 


1 +h- 


^ + HOT 


Thus condition (15) implies that the stretching is quasi— uniform. 

Turkel considers several finite difference and quadrilateral based finite volume 
formulations and his main conclusions are that with quasi-uniform stretchings all second order 
techniques retain their accuracy locally. If the stretchings are exponential or faster, central 
difference approximations will usually deteriorate to first order accuracy, and in the case of a 
cell centered approximation, the finite volume approach can yield an inconsistent scheme. 

Another illuminating investigation of truncation error has been provided by Roe [67]. 
He considers vertex based finite volume schemes and examines the error that arises from using 
a trapezoidal integration rule to evaluate the flux integral around different cell types. He 
concludes that only a few specific cell types retain second order accuracy locally. In 
particular, among quadrilaterals only a parallelogram can admit a locally second order accurate 
approximation. Furthermore he shows that a triangular mesh can only have first order 
accuracy. 

This rather alarming conclusion would appear to put unstructured meshes at a serious 
disadvantage. However, recent work by Giles [36] suggests that although the truncation error 
on an unstructured mesh is first order, the solution error remains second order. Further 
investigation is required to shed more light on this difficult area. It is clear, however, that 
future investigations of accuracy must attempt to follow Giles' lead in estimating the solution 
error and not basing conclusions on the order of the truncation error alone. 

5.2 Wave Propagation 

Consider the scalar wave equation 

u +cu =0, c>0 (16) 

l X 


which represents the propagation of disturbance traveling to the right with constant velocity c. 
We now discretize the spatial terms to obtain the semi-discrete form 
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( 17 ) 


ft - c ( v - v 

ar x j+1 - Xj — j 

On a uniform mesh Xj + j — Xj = h = const, and this equation reduces to 

du i . (u i+i - "i-P 

art = ~ c 2h 

Equation (16) admits traveling wave solutions of the form 

i(kx - tot) h 

c 


(18) 


If we look for similar solutions of the finite difference approximation (18) we obtain a 
dispersion relation 


co = | sin(kh), -n < kh < 7t 

In other words, the mesh acts like a dispersive medium so that a component with wave number 
k will travel with a phase velocity 


c* = c 


sin(kh) 

IcH 


and a group velocity 

= c cos(kh) 


We observe that for small values of the wave number k, the group velocity is close to the 
continuum solution velocity c. For higher wave numbers this is no longer the case. In fact, 
when kh = rc/2 corresponding to a wave length 



4h 


the component has zero group velocity, and at the highest wave number kh = 7t, which 
corresponds to the shortest wavelength 


X = 2h 


the component travels upstream at speed c. 

An investigation of the wave propagation characteristics for different time integration 
schemes has been pursued by Giles and Thompkins [37], Trefethen [91] and Vichnevetsky 
[93—95]. Of particular interest is recent work by Vichnevetsky [95] who considers a 
non-uniform mesh and looks for solutions of the form 

u-A^t 
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The wave number k and amplitude are now slowly varying functions of x. As the wave 

number changes on an expanding mesh, it is possible for the group velocity to decrease, reach 
zero and become negative (see figure 15). A solution component can therefore undergo a 
reflection and, if the mesh expands in all directions, it is possible for the wave to become 
trapped. 
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~aa"" ^ — 




X 


Fig. 15. Wave Reflection on an Expanding Mesh 
(taken from Vichnevetsky, Ref. 95) 


Vichnevetsky’s analysis and indeed the results of all the investigators in this area is 
restricted to non— dissipative approximations. When artificial dissipation is present, these 
waves are damped and no undesired transients are left trapped in the mesh. However, the 
wave propagation analysis could have important implications for unsteady flow problems 
where the artificial dissipation should be kept small and not interfere with the true time 
accurate behavior. 
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6. Adaptive Meshing 


It has long been recognized that a mesh determined independently of the flow solution 
is unlikely to resolve all the flow features. Specific regions of the flow field will be subject to 
rapid variations in flow variables. Solutions of the Euler equations require the accurate 
capture of shockwaves and the resolution of vortices that can emanate from sharp edges in 
three dimensional computations. With the advent of Navier Stokes calculations it becomes 
necessary to ensure a fine mesh in wake regions and areas of separated flow. It is clearly not 
possible to achieve the necessary mesh resolution a priori without introducing an unacceptably 
large number of points throughout the flow field. The use of an adaptive mesh that evolves 
with the flow solution is therefore an unavoidable requirement. In fact, it is anticipated that 
an adaptive meshing capability will be an integral part of all three dimensional flow codes 
within a few years. 

6.1 Mesh Redistribution 

An option that received early attention is use of mesh redistribution to achieve the 
required density of mesh points. One dimensional mesh stretchings that bunch the mesh lines 
have been considered, and stretchings in more than one direction have also received attention 

[29.39.58.90] . Within the framework of numerical generation methods, it is possible to link 
the source and repulsion terms to some measure of the truncation error or variation in flow 
variables. Other numerical procedures based on a variational approach have been tried 

[17.90] . Here some flow property is included in the functional to be minimized. The 
Euler— Lagrange equation that one obtains can then be solved numerically to generate the new 
mesh. A common difficulty of all the methods based on mesh redistribution is the need to 
ensure that the mesh cells do not become too distorted and, in particular, do not cross over 
each other. It is therefore necessary to monitor very closely changes in cell size and shape, 
and to halt the mesh adaptation if tolerances on these constraints are exceeded. 

6.2 Mesh Enrichment 


An alternative approach which does not suffer from these difficulties is mesh 
enrichment. Here the mesh is successively refined whenever the truncation error, or some 
convenient measure of the flow variables, becomes too large. Notable work in this area is that 
of Berger and Jameson [14] and Dannenhoffer and Baron [25,26], The main difficulty of this 
approach is the need to keep track of the interfaces where refinement occurs. This requirement 
adds an overhead to the data handling and complicates the flow solver. It is also necessary to 
ensure conservation and stability at the interfaces where refinement occurs. These problems 
are difficult to treat in two dimensions; in three dimensions the complexity of this task is 
increased considerably. 

It is in the context of mesh adaption that unstructured meshes enjoy a distinct 
advantage over structured meshes. A triangular mesh can easily be refined without the need to 
introduce interfaces [1,23,41,42,53,54,57,60,61,63,83]. In other words the mesh remains 
entirely transparent to the flow solver. The addition of mesh points can be carried out either 
by explicitly refining triangles, or alternatively by retriangulating the mesh to include the extra 
mesh points. This latter approach is particularly well suited to the Delaunay triangulation 
technique and has been demonstrated very effectively by Holmes et al. [42]. Figure 16 shows 
the mesh in a channel with a bump for a supersonic onset flow. A Delaunay triangulator was 
used with adaptive mesh enrichment, and the concentration of triangles which formed around 
the shockwaves is clearly evident. An example of the former approach has been provided very 
recently for a three dimensional triangulation by Periaire [63]. 
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Fig. 16. Mesh Enrichment on a Triangular Mesh 
(taken from Holmes et al.. Ref. 42) 


7. Whither Mesh Generation ? 

It is a reckless though irresistible prospect to speculate on the future development of 
mesh generation methods. This is particularly apposite since we are now at a stage where the 
primary goal of treating complete aircraft configurations has been achieved. What remains is 
to bring these ’primitive’ methods to a stage where it is possible to handle shapes of arbitrary 
generality with ease, without the need for special treatment, or undue interaction from the user 
whenever a different configuration arises. Of particular importance will be the need to achieve 
the most accurate solution possible for a given number of points and, if feasible, provide an 
estimate of the accuracy obtained. The development of efficient and robust adaptive meshing 
techniques is therefore a mandatory requirement. Not only should extra points be added to 
refine the mesh where needed, but it is also necessary to ensure that the shape and, in 
particular, aspect ratio of the cells remains acceptable. This may require the addition of points 
in some areas and the deletion of points in other regions. A deeper understanding of the 
dependence of the solution on mesh stretching and distortion is therefore required. Techniques 
that will either enrich or prune the mesh in selected regions must also be implemented as part 

of an adaptive mesh generation package. It is the author's opinion that it will be the 
unstructured meshes that best meet these stringent requirements for generality and adaptability. 

However, much effort has been devoted to the development of structured mesh generation 
schemes and it is unlikely that they will be abandoned in the near future. Indeed, the future 
vigor of research into mesh generation can best be served by the competitive pursuit of both 
approaches. 
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